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I.   INTRODUCTION 

Recent  research  at  the  NPS  in  atmospheric  global  predic- 
tion has  been  primarily  concerned  with  the  five-level  global 
primitive  equation  model  initially  designed  and  programmed 
by  Dr.  F.J.  Winninghoff  and  further  developed  by  Elias  (1973), 
Mihok  (1974),  McCollough  (1974)  and  Maher  (1974).   Their 
numerical  experiments  proved  valuable  but  pointed  out  the 
inflexibility  of  the  model  for  research  purposes. 

Arakawa  and  Mintz  (1974)  described  an  atmospheric  global 
prediction  model  with  a  new  finite  difference  formulation 
and  greater  vertical  resolution.   The  model  employs  a  new 
horizontal  distribution  of  variables  to  improve  the  geo- 
strophic  adjustment  process.   The  purpose  of  this  research 
was  to  program  and  check  out  a  modified  version  of  the  new 
Arakawa  scheme.   The  fundamental  approach  was  to  develop  a 
flexible  program  for  use  in  a  number  of  projects  including 
modeling  of  tropical  circulation.   Eventually  the  model's 
performance  can  be  compared  with  the  performance  of  FNWC's 
global  model. 

The  model  was  initially  developed  with  two  levels  and  a 
variable  grid  to  aid  in  debugging.   The  atmosphere  was 
treated  as  adiabatic  and  frictionless  and  the  water  vapor 
continuity  equation  was  not  considered.   The  initial  condi- 
tions for  all  experiments  were  analytically  derived  in  the 
same  manner  as  Maher  (1974).   The  advantages  of  analytic 
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initial  conditions  were  significant  in  reducing  computer  time 
during  check  out  of  the  model. 
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II.   MODEL  DESCRIPTION 

The  differential  equations  are  essentially  the  same  as 
those  described  by  Arakawa  (1972).   The  integration  scheme 
was  carried  out  on  a  staggered,  spherical,  sigma  coordinate 
system.   For  an  adiabatic,  frictionless  model  the  primitive 
equations  form  a  closed  set.   No  sink  or  source  terms  were 
considered  in  the  course  of  development  or  check  out. 

A.   VERTICAL  COORDINATE  SYSTEM 

The  model  uses  the  non-dimensional  sigma  coordinate  sys- 
tem as  described  in  Haltiner  (1971).   The  two  levels  divide 
the  troposphere  in  half  and  the  tropopause  is  assumed  iso- 
baric.   The  sigma  coordinate  is  defined  as 

P  -  Pt 
o   =        -  Z    ,  (2.1) 

where  p  is  the  pressure,  p.  the  constant  tropopause  pressure 
and  it  is  the  terrain  pressure.  The  terrain  pressure  is  fur- 
ther defined  as 

tt  E  ps  -  pt  ,  (2.2) 

where  p   is  the  surface  pressure.   It  follows  from  equation 
(2.1)  that 


a  =  0  at  p  =  p   , 


a  =  1  at  p  =  ps  , 


(2.3) 


which  are  the  vertical  boundaries  of  the  coordinate  system 


16 


(7=0    (7=0 


(7=.2  5 


>> 


*%y 


(7  =  .5 


(7  =  .75 


(7  =  1     (7=0 


Figure  1.   The  sigma  (a)  coordinate  system  as  used  in  the 
model.   p.  is  200  mb. 
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The  boundary  condition  at  a=l  is  a=0  and  at  a=0  it  is  as- 
suraed  that  a=0  where  a  is  defined  as  da/dt. 

B.   PRIMITIVE  EQUATIONS 

The  continuous  form  of  the  primitive  equations  is  de- 
scribed in  the  orthogonal  curvilinear  coordinates  £  and  n 
where  £  =  A  (longitude)  and  n  =  <J>  (latitude).   An  area  ele- 
ment is  — A£Ari  and  the  lengths  of  the  sides  are  —  and  —  , 
mn  &  m      n 

where 


—  =  a  cos  d>  and  —  =  a. 
m         r     n 


The  difference  element  is  illustrated  below, 


rj  +Arj 


£+A$ 


The  mass   continuity   equation    is 


8    ,  ir  ,.         8     ,    u,         9    f   v.         8    .irav    _    n 


(2.4) 
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The  ^-component  of  the  horizontal  equation  of  motion  is 

3  ,    TT   N     3   ,TTU  N     3  -,1V  N     3  ,TTO      . 
rr( U)  +  -7r-p( U)  +  -r— ( U)  +  -r— ( U) 

3tvmn  '    3£v  n      3nv  m  y    3avmn  ' 

r  f      ,31       3    1N,        TTr3$  3lT,      IT   -p        /o  cn 

-  [ —  +  (v-Trp-  —  -  u^ )]ttv  +  —  [tt=  +  aa  t-t]  =  —  F_  .   (2.5) 

mn    v  3£  n     3n  m       n  3?      3£    mn   £ 

The  n-component  of  the  horizontal  equation  of  motion  is 

3  ,    TT   .     3   ,TTU  v     3   /TTV  v     3   .TTO   N 

It(S5v)  +  Tf<17v)  +  ¥H(1TV)  +  I?(S5V) 

t  t J.  +  (vlfi-  »|-i)]wu  +5[|4  +  oa|S]  -J-F   .   (2.6) 

mn     3£  n     3n.  m       m  9n      3n    mn   n 


The  first  law  of  thermodynamics  is  written  in  the  following 
manner : 

S  i  =  uo  +  Q-  (2-7) 


where 


00  =  ft  =  ™  +   a(ft  +  V'V)tt  (2.8) 


and   the   flux   form   is 

-^(TrCpT)    +    Vq-(ttV   CpT)    +  |_(TTaCpT)    =   Tr(ooa    +   Q).  (2.9) 


Using  the  relationship  for  potential  temperature, 


T  =  6(£-)K  ,  (2.10) 

Po 


the  final  form  in  curvilinear  coordinates  is 
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3    (JL  c   T)    +  i^Si  C   T)    +  f-(^  C  T)  +    (2-)*  i-(™   C    9) 

3t   mn      p    \         Hv  n       p             3nv  m      p    '  vp0         aavmn      p    J 

=   naa  ttC— -)    +  -  -57-  +  -  3-      +  —  Q    .  (2.11) 
dt   mn          n    dt,        m   dri          mn 


The  equation  of  state  is 

a  =  —  (2.12) 

P 

and  the  hydrostatic  relationship  is 

6$  =  -TTct6a  .  (2.13) 

A  complete  set  of  symbols  for  the  above  equation  may  be 
found  in  the  front  of  this  report. 

C.   VERTICAL  INDEX 

The  index  k  is  used  to  identify  the  vertical  levels.   At 
the  upper  boundary  k=0,  p=p.  and  at  the  lower  boundary 

k=K+l  and  p=p  .   The  variables  V  and  T  are  carried  at  the 

*      s 

odd  levels  while  tto  is  carried  at  the  even  levels.   The  fol- 
lowing definitions  are  required  to  complete  the  vertical 
system: 

Aak  ~   ak+l  "  ak-l  '  (2'14) 

K 
^  Aak  E  1  ,  (2.15) 

where  the  summation  is  over  odd  k   (see  figure  2). 
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INDEX  COMPUTED  SIGMA 

k  VARIABLES 

0 77"  (J— Q       <j  =  0 

1 VTg" 

_2 7T_&__  <J-  o2 

_3 VT# 


lc-2 VT<2>" 

k-1 ' rr  a  <7—  a.   4 


VT^ 


k-H 7T  (7  <7  =  <7.  +  1 

k+2  vjJIL 


K-2 

VT<? 

K-1 

7T<7 

K 

VT  #" 

K+1 

7T  a  — 

o  =  cK_, 


VERTICAL   INDEX 

Figure  2.   In  the  above  figure  k  is  a  variable  vertical 
index,  sigma  (a)  is  the  dimensionless  vertical  coordinate, 
V  is  the  horizontal  vector  velocity,  tt  is  the  terrain  pres- 
sure,  a  is  the  vertical  velocity,  T  is  the  temperature  and 
$  is  the  geopotential. 
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D.   HORIZONTAL  GRID  AND  DISTRIBUTION  OF  VARIABLES 

The  selection  of  the  horizontal  distribution  of  variables 
is  a  function  of  two  distinct  processes.   The  first  is  proper 
simulation  of  the  geostrophic  adjustment  process  and  the 
second  is  proper  simulation  of  slowly  changing  quasi-geo- 
strophic  motion  after  it  has  been  established  by  geostrophic 
adjustment.   As  shown  by  Winninghoff  (1968)  geostrophic  ad- 
justment depends  on  how  the  variables  are  distributed  over 
the  grid  points.   Winninghoff  used  the  following  equations, 
which  are  the  simplest  ones  in  which  geostrophic  adjustment 
can  take  place,  to  demonstrate  five  possibilities  for  the 
placement  of  the  dependent  variables: 

f--fv   +   gf  =   0,  (2.16) 


g   ♦   ?u    ♦   gff   -    0    ,  (2.17) 


8  ♦  *<1i  ♦  H>  -  °  •  <2-18> 


The  above  equations  represent  an  incompressible,  homogenous, 
non-viscous,  hydrostatic,  rotational  fluid  with  a  flat  bot- 
tom and  a  free  surface,  where  u  and  v  are  the  velocity  com- 
ponents, h  is  the  depth  of  the  fluid,  f  is  a  constant 
Coriolis  parameter,  t  is  the  time,  x  and  y  are  the  horizon- 
tal coordinates  and  g  is  gravity.   The  five  possible  distri- 
butions of  the  dependent  variables  h,  v  and  u  are  shown  in 
figure  3.   Scheme  B  was  used  in  the  Mintz-Arakawa  two-level 
global  circulation  model  (Langlois  and  Kwok ,  1969)  and 
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Figure  3.   The  above  figure  represents  the  placement  of  the 
dependent  variables  h,  v  and  u  where  h  is  the  depth  of  the 
fluid,  u  and  v  are  the  velocity  components  and  d  is  the  grid 
distance  (for  the  simplest  case  of  geostrophic  adjustment). 


23 


scheme  E  is  used  in  the  global  model  now  under  further  de- 
velopment at  FNWC  (Maher,  1974). 

In  the  one-dimensional  case  it  was  shown  that  both 
scheme  B  and  C  adequately  simulated  geostrophic  adjustment. 
In  the  two-dimensional  case,  however,  scheme  C  was  shown 
most  satisfactory  to  simulate  geostrophic  adjustment,  except 
where  A  /d  is  less  than  or  close  to  one.   The  quantity  X0  is 
the  Rossby  radius  of  deformation  and  d  is  the  grid  distance 
(Arakawa  and  Mintz,  1974).   The  Rossby  radius  of  deformation 
is 

xR  =  /gT/f 

where  H  is  the  mean  value  of  h,  g  is  the  acceleration  of 
gravity  and  f  is  the  constant  Coriolis  parameter.   The  con- 
dition XR/d  <  1  is  an  abnormal  case,  therefore,  the  horizon- 
tal distribution  of  variables  is  based  on  scheme  C.   (See 
figure  4. ) 

E.   TIME  DIFFERENCING  • 

The  time  differencing  is  carried  out  in  thirty  minute 
sequences.   The  initial  step  in  each  sequence  is  a  two  part 
Matsuno  scheme  represented  by  the  following  notation: 

F*  =  Ft  +  At  ~   (Forward),  (2.19) 

o  t 

Ft+1  =  pt  +  At  ||  (Backward)  .  (2.20) 

o  t 
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Figure  4.   The  above  figure  represents  the  horizontal  distri- 
bution of  dependent  variables  with  polar  modification.   P  in 
the  ij  index  system  represents  the  north  pole,  1  represents 
the  south  pole  and  i  represents  a  meridian.   tt  represents 
the  variables  T,  <!>  and  tt  carried  at  "TT-points."   u  and  v  are 
the  horizontal  components  of  velocity. 
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F  is  a  vector  representing  the  dependent  variables.   The 
superscript  t  represents  the  time  step,  the  superscript  * 
represents  the  results  of  an  intermediate  step  and  At  is  the 
time  interval  of  a  single  step.   The  remaining  time  steps  in 
a  sequence  use  the  leapfrog  scheme  as  follows: 

Ft+1  =  Ft_1  +  2At  ||   .  (2.21) 

At  the  end  of  each  sequence  provision  is  made  for  the  calcu- 
lation of  the  source  terms.   The  present  model  doesn't  in- 
corporate source  term  calculations. 

F.   AVERGING  THE  PRESSURE  GRADIENT  AND  ZONAL  MASS  FLUX  NEAR 

THE  POLES 

A  problem  that  may  arise  at  higher  latitudes  is  computa- 
tional instability.   Computational  instability  is  a  result 
of  the  convergence  of  the  meridians  to  the  poles  causing 
greatly  reduced  grid  distances  along  a  latitude  circle. 
Therefore,  some  technique  must  be  employed  to  eliminate  the 
stability  problem.   One  technique  involves  reducing  At  in 
higher  latitudes  but  it  would  impose  serious  programming 
difficulties  and  require  more  computer  time.   Another  tech- 
nique would  be  to  reduce  the  number  of  grid  points  in  a 
latitude  circle  at  higher  latitudes,  however,  it  also  would 
impose  programming  difficulties.   A  third  technique  is  zonal 
smoothing  which  was  developed  by  Arakawa  and  Mintz  (1974). 
Smoothing  preserves  the  integrity  of  the  horizontal  grid 
while  maintaining  a  constant  time  step.   Zonal  smoothing  is 
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used  in  this  model  and  also  by  FNWC  to  eliminate  computa- 
tional instability  at  high  latitudes. 

Zonal  smoothing  involves  expanding  the  pressure  gradient 
and  zonal  mass  flux  into  a  Fourier  series  and  reducing  the 
amplitude  of  each  wave  component  by  a  factor 

q  =    D  cos  $  (<>    99s 

b   C0At  sin  (rod)  {4.4* ) 

where 

S  =  Stability  Coefficient, 

D  =  Grid  distance  at  the  equator, 

Co  =  Phase  speed  of  the  fastest  gravity  wave, 

At  =  Time  step, 

m,  =  Wave  number , 

d  =  Grid  distance  in  degrees. 

If  the  value  of  S  is  greater  than  one,  smoothing  is  not  re- 
quired (Arakawa  and  Mintz,  1974).   As  pointed  out  by  Arakawa 
the  smoothing  operation  does. not  smooth  the  fields  of  vari- 
ables, because  it  is  simply  a  generation  of  multiple  point 
difference  quotients.   The  bar  operators  in  Chapter  III  in- 
dicate zonal  smoothing. 

G.   PROGRAM  FORMAT 

The  program  format  was  patterned  after  the  Mintz-Arakawa 
two-level  global  model  as  described  by  Gates,  et  al.  (1971) 
and  Langlois  and  Kwok  (1969)  and  modified  to  include  scheme 
C  distribution  of  variables  and  a  Matsuno-leapf rog  time  in- 
tegration scheme.   Figure  5  gives  a  simplified  flow  diagram 
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Figure  5.   A  simplified  flow  diagram  for  5-day  forecast 
using  analytic  initial  conditions. 
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for  a  five-day  forecast.   The  main  program  controls  the 
overall  time  period  of  the  forecast  as  well  as  input  and 
output.   Time  integration  for  each  30-minute  sequence  is 
controlled  by  subroutine  Timestep.   The  forcing  functions 
are  calculated  in  Comp  1  and  Comp  2  with  zonal  smoothing  ap- 
plied as  necessary.   The  remaining  subroutines  are  peripheral 
to  the  main  flow  of  the  model. 

A  main  consideration  in  the  program  was  flexibility, 
therefore,  the  vertical  structure  along  with  the  horizontal 
grid  are  variable.   The  model  can  be  used  to  simulate  sever- 
al different  horizontal  and  vertical  structures  including 
the  FNWC's  5- level  global  model. 
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III.   FINITE  DIFFERENCE  EQUATION 

The  finite  difference  equations  given  in  this  chapter 
were  developed  by  Arakawa  and  Mintz  (1974).   The  advection 
terms  of  the  equation  of  motion  were  modified  to  eliminate 
the  diagonal  flux  calculations.   The  short-term  performance 
should  not  be  affected  by  elimination  of  the  diagonal  flux 
terms  and  considerable  computer  time  was  saved. 

The  notation  used  in  this  chapter  represents  the  specif- 
ic variable  centered  in  the  ij  index  system.   Figure  6  con- 
tains an  example  of  "iT-centered"  and  "u-centered"  notation. 
Ti-centered  notation  is  used  for  the  continuity  equation  and 
the  thermodynamic  equation  in  which  the  dependent  variable 
is  carried  at  a  7T-point.   U-centered  notation  is  used  for 
the  ^-component  of  the  equation  of  motion  and  v-centered 
notation  is  used  for  the  n-component  of  the  equation  of 
motion.   Refer  to  the  front  of  the  report  for  a  complete 
list  of  symbols. 

A.   CONTINUITY  EQUATION 

The  following  form  is  used  for  the  continuity  equation 
given  in  equation  (2.4): 


*'j  +  Fk    .  -  Fk  ,  .  +  Gk  ...  -  Gk  .  1 

at      i+i,j    k-£,j    i,j+i    i,j-i 

+  1  (sk+1  -  sk-x)  =  o  ,  (3.1) 

Aak   i,j    i,j' 
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^CENTERED 
NOTAT  ION 


Figure  6.   An  example  of  the  notation  used  to  describe  the  finite  dif- 
ference equations.   The  continuity  equation  is  described  on  a  "Tr-centered" 
grid  in  which  the  F  and  G  symbols  are  flux  calculations  in  their  re- 
spective directions.   The  "u-centered"  grid  is  an  example  used  to  de- 
scribe the  ^-component  of  the  equation  of  motion  where  F      and  g  are 
also  flux  calculations. 
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where 


F 


U.i   =  ""T^iH.J  ('i+i,j  +  "i.J>  '  <3"2> 


(bar  indicates  zonal  smoothing) 


i.J+i  "  *(Vm  >1-J+|  (*i-J+l  *  *i§j)  ■  <3'3> 


n.   =  ¥.  .  (MM) 

i,J     i,J    mn  '   . 

-1-  >  J 


S.    =  n.    a. 

1,J    1,J  1,J 


The  tendency  equation  is 


at     k=i  ui+i,j   *i-i,j  '  ijH   i,j+£'    •  v   ; 


The  vertical  motion  equation  is 
k+1     l«  /T1k       ^k     ,  „k       ^,k    NA  k 


S 


i,j    k=l   1+4, j    l-i, j    1,3+i    i,J-4 


k+1  ^U.  (3.5) 

at 
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B.   PRESSURE  GRADIENT  TBftM  AND  HYDROSTATIC  EQUATION 
1.   Hydrostatic  Equation 

The  geopotential  for  a  level  K  is  given  by 

K  RT 

*K   _*S       r  '  k    i.i.k 

$  .    =  $  .    +  ,  ^^  it  .  .a  — r-^  Act 

i,J    i,J   k=l  i,j     k 


K-2,   .  ,1     „,...    _k+2  k       k  k 
k  P   1.0    Po    1;j     P0  ij 

g 

where  $   is  the  geopotential  at  the  surface, 

~k+l    £n  6k  -  £n  6k+2 

0     =  -j  j  . 


Qk+2    Qk 

Given  $  ,  the  geopotential  for  the  remaining  levels  is  found 
using  the  following  relationship: 


$k  .  -  $k+2  =  C  [(£ )     -  (£-)    ]  6k+1  .       (3.7) 

i.J     iJ     P   Po  \fj  Po  ij    X>J 


2 .   n-component  of  the  Pressure  Gradient  Force 

The  pressure  gradient  force  in  the  meridional  equa- 
tion of  motion  as  given  in  equation  (2.6)  is 

"  m1  IK*!  i+Z    +  *1  i  *><*?  i+*  "  *i  i  *> 

m-j  l.J  +  2        1,J~2      l,J  +  2        1,J~2 

+  ((Traa)k  ,  ,  +  (Traa)k    ,)(ir.  .  ,  -  tt   .  ,)]  ,    (3.8) 

l,J  +  2  l.j~2      l-J  +  2        l.J~2 
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where 


k  k  RTi  i+x 

(Tract).  .,,  =  it.  ...  a  ^ — a 

'i,j+£    i,j+i     p,  i+A 


3.   ^-component  of  the  Pressure  Gradient  Force 

The  pressure  gradient  force  in  the  zonal  equation  of 
motion  as  given  in  equation  (2.5)  is 


M     itU.^   •  +  7T.   .)<•*.   ■  "  4>k   .) 

n.      i+l, J     i.J    i+l, J     i,J 

•J 


+  ((Traa)k    •  +  (Traa)k  .)(tt..-  .  -  it  .  .)]   .       (3.9) 
'i+l,j         i,J    i+l, J     i,J 


The  zonal  momentum  flux  as  given  in  equation  (2.5)  is 


-Fi-i,j(ui,J +  ui-l,j)k  +  gi,j+i(ui,j+l +  ui,/ 


u    ,                              xk,  1   lrMi,k+l,  k+2  ,   k   >. 

K.  .  -,  (u .  .+u.    -,  )  1  +  — r  4 1 S  .  .   (u.    +  u .  .) 

Bi,j-£k  i,j     i,J-l'  A„k    i,j     i,j     i,j' 

,u,k-l,  k      k-2N, 
!.  '  .   (u .    +  u.  . )]  , 
i,J    i,J    i,J  J 


Aa 
■  ^'k-1(nk        +  uk"2)l  ,  (3.10) 


where 


ni  i  E  i(ni+i  i  +  ni-i  i>  ' 

1  »  J  1^2  ,  J        1   2  ,  J 


34 


,u 


S   .  E  *(S.^,  .  +  S.  ,   .) 


F  and  g   are  defined  as  follows 


FU,3   E  i(F*iH,J+i  +  **U.i  +   Ft+i,j-i>  •         (3-11) 


U  *  *  *  * 

..,  =    iCG-.j.  .  +  G.,,  .,-  +  G.  ,  .  +  G.  ,  ....)  ,   (3.12) 


where 


F*  .  e  £(F.,,  •  +  F.  ,  .)  ,  (3.13) 

i,J       i+i,J     i-i,J 


G*  .  E  i(G.  ..,  +  G.  .  x)  .  (3.14) 

i,J       i,J+i     i,J-i 


The  meridional  momentum  flux  is  similar  to  equation  (3.10) 
with  u  replaced  by  v.  F   ,  g  ,  II   and  S   take  the  following 
form: 


FV^x     ■    =   iCF.j.-,  -..I  +  F.  ..,  +  F.  .  i  +  F.,..  .  ,)  ,   (3.15) 


■I.J+4 '-  i(Gi+l,j+i +  2Gi..j+i +  Gi-1,H» .  <3-16) 

nI.J  E  i(ni,j+i  +  "i.J-i>  •  (3-17) 

§I,J E  i(§i,M +  *l.J-»>  •  (3-18) 


F   is  defined  by  equation  (3.13)  and  G   by  equation  (3.14). 
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D.       CORIOLIS   FORCE 

A  variable  C.     .is   defined   at   Tr-points   by 


ck  .  e  f .  .(AMll)  _  i(u    .  +  u    ,}k 


[(f)     -(f)    1  , 
m  j+£     m  j-i 


where  f .  .is  the  Coriolis  parameter. 
i  >  J 

The  zonal-component  of  the  Coriolis  force  as  given  in 
equation  (2.5)  is 


[*, 


k      ,  ,k 


*1+*.J  Ci+i,J  (Vi+i,J+i  +  Vi+i,j-i) 

V*  i  Ci_i  i  (vi-^  i+*  +  vi-i  i-^)k]  '         (3-19) 

J-   2  >  J     J-   2  >  J       -L   2>JT2         J-   2  >  J   2 


The  meridional-component  as  given  in  equation  (2.6)  is 


-i[ir     ,  C   .  ,  (u.  x    .  i  +  u._4  -i) 

-»-  >  J   2     1>J   2       1T2  )  J   2         J-   2  ,  J   2 


+  ir.  ..,  C*      (u  '  ,    ,  +  u  x        j,)1"]  .  (3.20) 

l,J+2     l,J+2       l+2,J+2         1~2, J+2 


E.   THERMODYNAMIC  ENERGY  EQUATION 

The  thermodynamic  equation  corresponding  to  equation 
(2.11)  may  be  written  as 
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|r  <n.         Tk    .)    +  Fk  (Ti+lJ    !  Tir,1)k 

-   Fk_,         (Ti>J    ^Ti-1-1)k   +   Gk  (Tl>J+12+  Ti>J)k 

1       2  ,  J  &  1  ,  J +  2  ^ 

_Gk  (Ti,j  *Ti,.1-i)k  +  ^-[sk+1  (^i)K  ek+1 

•k-1       Pi      i    K    ~k-1  1  k  ^i     i 

-Si,j    <p7^>      9i,j=    =^'^)i,j3T^ 


1T2  ,  J 


+    (u^r)  ((7raa)k    .    +    (-rraa)k  .     .)(tt.     .    -   tt  .    .     .) 

n  %    i    -i  i,J  'i-l,j/v    i,j  l-l, j' 


1-4,  J 

k 


+    (v^)  ((iraa)"    .    -    +   (iraa)?      )(tt  -   it        ) 

^    H+i  i>J+1  i,j        i,J+i  i|j 

1  >  J^2 


m 


+  (v--f)      ((Tract).  .  +  (Tract).  .1)(tt.  .  -  tt  .  .  .,  ) } 
™  s  -i  i       iJ        'i,j-l'v  i,j     i,J-ly 

1  >  J~2 


+  n.  .  Qk  .]  .  (3.21) 


F.   POLAR  MODIFICATION 

The  poles  of  a  spherical  coordinate  system  are  singular 
points  and  polar  velocity  components  cannot  be  defined. 
Terrain  pressure  (tt)  at  the  poles  can  change  as  a  result  of 
meridional  flux  at  all  points  P-i  and  P+i  where  v  is  carried 
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The  continuity  equation  is  modified  at  the  poles  by  omitting 
all  undefined  flux  terms.   Each  pole  is  treated  as  a  series 
of  points  in  order  to  facilitate  programming  (see  figure  7). 
Each  index  i,P  represents  the  shaded  area  shown  in  figure  8. 
The  continuity  equation  is  integrated  at  each  point  repre- 
senting the  pole  and  then  averaged  to  determine  a  polar  val- 
ue.  The  thermodynamic  equation  and  the  vertical  velocity 
are  treated  in  a  similar  manner.   The  advective  terms  in  the 
equation  of  motion  are  given  special  treatment.   The  remain- 
ing terms  in  the  equation  of  motion  are  unchanged  except  in 
the  case  where  they  are  undefined  and  omitted. 

1.   The  Advective  Terms  in  the  ^-component  of  the  Equa- 
tion of  Motion 

The  polar  modification  of  the  ^-component  of  the 

equation  of  motion  is  shown  in  figure  9.   It  was  modified 

from  the  Arakawa  scheme  by  eliminating  the  diagonal  flux 

terms.   The  following  form  is  used: 

9(nuuk)  +  if  P*u  ,  .k 

at        i,p_i     2[Fi+i,P-i  (ui,p-i      ui+i,p) 

r-*U  /  Nk  u  ,  vk, 

-     ^i-i,P-l    (ui,P-l   +   ui-l,P-l>      -    Bi.P-Va    (Ui,P-l   +   Ui,P-2}    ] 

1      lrMi,k+l,    k+2  k  N        Mi,k-1,    k  k-2      v, 

+  7k  **§        (ui,p-i  +  ui,p-i)  -  s        (ui,p-i  +  ui,p-i)]  ■ 

Aa  '  . 

(3.22) 

where 


*u  i  *  * 

Fi-i,P-l   =   4(4Fi-ir,P-l  +  Fi-i,P-2}    ' 
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Figure  8.   Each  index  i,P  in  the  continuity  equation  is 

represented  by  the  shaded  area.   it  at  the  poles 
can  only  change  as  a  result  of  G.   The  thermo- 
dynamic equation  and  vertical  velocity  are  treated 
in  the  same  manner. 
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Figure  9.   The  polar  modification  of  the  u  equation  of  motion 
F       and  g   are  flux  terms.   The  shaded  portion 
represents  the  area  associated  with  each  variable 
u. 
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g   given  by  equation  (3.12) 


rU 


ni,P-l  '    *<ni-4,P-l  +  ni+i,P-l> 


u 


s: 


i,p-i  ~   2^si-i,p-i  +  si+i,P-i} 


2.   The  Advective  Terms  in  the  n-component  of  the  Equa- 
tion of  Motion 

The  polar  modification  of  the  n-component  is  repre- 
sented in  figure  10.   It  was  modified  from  the  Arakawa 
scheme  in  a  similar  manner  as  equation  (3.22).   The  follow- 
ing form  is  used: 

9(nvvk)     +  i,  F*v    ,  +      ,k 

at     •  p  A  +  2l  h i+i,P-i(Vi+l,P-i   Vi,P-i} 

-L  ,  r—  2 

-  /C*,P-*(vi,P-i  +  vi-i,P-i)k  -  gI,p-i(vi,P-i  +  Vi,p-V2)k] 

+  _JL_  irQV,k+l   k+2      k     v    §v,k-l   k        k-2   ., 

+     V  2lb        ^V-i  P-i     Vi  P_A;  "  b        ^Vi  P--1-     Vi  P-i;J   » 

«   K  1,^-2         l,-r_2  1,-t^   2         1,-t^   2 

(3.23) 
where 

P*v      =  £(F*       +  F*     ) 
^±-h,P-h    ~    *K    i-l,P-l    i,P-l;' 

v 
g   given  by  equation  (3.16), 

nI,P-i  E  "(IIi,P  +  ni,P-l}  ' 


SijP_i  =  2(S.p  +  s.p_1)  . 
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Figure  10 


The  polar  modification  of  the  v  equation  of 

r-.  *  V  *  V 

motion.  F  and  g  are  flux  terms.  The  shaded 
portion  represents  the  area  associated  with  each 
variable  v. 
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IV.   ANALYTIC  INITIALIZATION 

The  initial  conditions  were  similar  to  those  used  by 
Maher  (1974).   Analytic  initial  conditions  offer  several  ad- 
vantages in  the  initial  evaluation  of  a  model;  they  simplify 
the  otherwise  difficult  task  of  balancing  and  interpolating 
initial  conditions  from  constant  pressure  surfaces  to  sigma 
(a)  surfaces,  they  allow  the  angular  phase  speed  to  be  esti- 
mated from  a  non-divergent  model  and  they  allow  the  simula- 
tion of  atmospheric  states  that  otherwise  would  be  difficult 
to  simulate. 

A.   ANALYTIC  BALANCING 

The  initial  terrain  pressure  and  velocity  fields  were 
obtained  from  the  stream  function  solution  to  the  linearized, 
non-divergent  vorticity  equation  (Haurwitz,  1940).   The  com- 
plete vorticity  equation  was  later  solved  by  Neamtan  (1946). 
The  stream  function,  \\> .    may  be  written  as  follows: 


mk 
\\)   =   A   sin(m,X   -   vt)    sin   <f>    cos        <j>   -   Ba2    sin  <J> ,  (4.1) 


where  A  and  B  are  constants,  v  is  the  angular  velocity,  m, 
is  the  wave  number  and  a  is  the  radius  of  the  earth.   The 
angular  phase  speed  is  given  by 

_v_  =  r  N(N+1)  -  2   _    2ft  ,4  2, 

rak       N(N+1)       N(N+1)  '  y    '     J 

where  v/m,  is  the  angular  phase  speed,  N=m,+1  and  ft  is  the 
angular  velocity  of  the  earth.   Harmonic  waves  defined  by 
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the  stream  function  will  move  with  a  constant,  angular  veloc- 
ity without  changing  shape  assuming  a  baratropic,  non-diver- 
gent atmosphere  (Haurwitz,  1940).   However,  the  equations 
used  in  this  model  are  for  a  divergent  atmosphere  and  will 
give  a  smaller  contribution  from  the  second  term  in  the  equa- 
tion (4.2),  especially  with  small  wave  numbers. 

Equation  (4.1)  is  used  as  the  forcing  function  to  obtain 
the  geopotential  from  the  non-linear  balance  equation 
(Phillips,  1959).   The  geopotential  perturbation,  $  ,  may  be 
written 

$'  =  a2A((J>)  +  a2B(cj))  sin  mkX  +  a2C(4>)(2  sin2mkA  -  1)   (4.3) 

where 

r  A  2    2mk 

^  a 


2m 


[(mk+l)  cos2<J>  +  (2mk  -  mR  -  2)  -  -~r^]     ,       (4.4) 


2(fi+B)A       m 
B<*>  ■  (m,+l)(m,+  2)  cos    *  [<mk  +  2mk  +  2) 


-  (m  +1)2  cos2  (J)]   ,  (4.5) 


a  2     2m, 
C(4>)  =  K-4")   cos   K  $  [(mk+l)  cos2c{>  -  (mk+2)]  .       (4.6) 
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B.   ANALYTIC  WINDS 

The  wind  components  were  obtained  as  follows: 

u     a  8<J>  »  ^*  '; 


v  =  ^ -^  M  8) 

a  cos  <}>  9X  *  v^-o; 


Performing  the  operations  indicated  above  the  u  and  v  com- 
ponents may  be  written 

i  V1 

u  =   -  —[A   sin(m.  X   -   vt)    cos  <f>  —  m,  A 

a  k  k 


m,  -1 


sin(m   A   -    vt)    cos  <j>    sin2cf>   -   Ba2    cos   <j>]     ,  (4.9) 


1  mk_1 

a  k  k 


C .   TEMPERATURE 

The  temperature  field  was  derived  in  accordance  with  the 
NACA  standard  atmosphere  as  follows  (Haltiner  and  Martin, 
1957): 

T(°K)  =  288  -  0.0065Z,  (4.11) 

Z  (meters)  =  44308  [1  -  (10^325)        3   ,   (4.12) 
where   Z  <  10769  meters. 
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V.   MODEL  PERFORMANCE 

The  evaluation  of  the  model  was  a  two-step  process.   The 
first  step  involved  the  use  of  a  coarse  grid  and  was  essen- 
tially a  scheme  to  minimize  computer  time  while  debugging. 
The  second  step  involved  the  use  of  a  finer  grid  for  better 
resolution.   The  majority  of  the  experiments  presented  in 
this  chapter  were  run  with  the  fine  grid. 

The  coarse  grid  consisted  of  16  points  N-S  which  gave  a 
12  degree  separation  between  grid  points.   Since  strictly 
analytic  initial  conditions  were  used  it  was  only  necessary 
to  integrate  over  one  wave  using  cyclic  continuity  and  there- 
fore a  10-point  E-W  grid  was  adequate.   The  E-W  grid  dis- 
tance was  a  function  of  wave  number.   Wave  numbers  4,  8  and 
12  were  used  in  the  experiments.  Wave  number  4  gave  an  E-W 
grid  spacing  of  9  degrees,  wave  number  8  gave  a  grid  spacing 
of  4.5  degrees  and  wave  number  12  gave  a  grid  spacing  of  3 
degrees.   The  fine  grid  used  46-points  N-S  which  gave  a  grid 
spacing  of  4  degrees.   The  E-W  grid  was  similar  to  the 
coarse  grid. 

In  all  cases  the  model  was  run  with  two  levels,  a  flat 
earth  and  no  source  or  sink  terms.   The  time  step  (At)  was 
six  minutes.   Fourier  analysis  of  the  surface  pressure  field 
was  used  to  compute  phase  speed  and  wave  amplitude  in  all 
but  experiment  III  in  which  a  Fourier  analysis  of  the  wind 
field  was  computed.   Calcomp  charts  were  plotted  at  12  hour 
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intervals  and  consisted  of  one  wave  of  the  surface  pressure 
field  over  the  complete  N-S  grid.   The  phase  speed  and  wave 
amplitude  were  initialized  through  the  constants  A  and  B 
given  in  equation  (4.1),  where  A  is  the  amplitude  of  the  dis- 
turbance and  B  is  the  amplitude  of  the  mean  flow. 

Experiment  I.   Wave  numbers  4,  8  and  12  were  used  to  de- 
termine the  phase  speed  over  the  coarse,  16-point  grid.   The 
phase  speed  was  set  at  10  degrees  per  day  by  adjusting  the 
amplitude  of  the  mean  flow  (B).   The  graphs  of  phase  angle 
(degrees  longitude)  vs  latitude  are  shown  in  figures  11,  12 
and  13.   The  observed  phase  speed  approached  20  degrees  over 
48  hours  for  all  three  cases. 

Experiment  II.   This  and  all  subsequent  experiments  were 
run  over  the  46— point  N— S  grid.   Extern?.!  gravity  waves  were 
simulated  with  wave  number  4  by  setting  the  mean  flow,  the 
u  and  v  components  of  the  wind  and  the  Coriolis  terms  to 
zero.   The  theoretical  period  of  the  gravity  wave  was  calcu- 
lated using  the  following  approximation: 

_  2(3.14159)  ,5  1, 

v  P 

where 


m    /N(N+1) 
p     a 


C  is  the  zonal  phase  speed  (set  at  300  m/sec),  K   is  the 
two-dimensional  wave  number,  N=m  +1,  a  is  the  radius  of  the 
earth  and  T   is  the  period.   The  above  approximation  gives  a 
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period  on  the  order  of  7  hours  for  wave  number  4.   The  ob- 
served oscillations  of  the  surface  pressure  field  as  a  func- 
tion of  time  are  given  in  figure  14.   The  period  is  on  the 
order  of  6-7  hours. 

Experiment  III.   In  this  experiment  pure  advection  was 
simulated  with  wave  number  4  by  removing  the  Coriolis  terms, 
the  pressure  gradient  terms  and  the  vertical  advection 
terms.   The  surface  pressure  and  the  temperature  were  kept 
constant.   The  u  and  v  components  of  the  wind  were  Fourier 
analyzed  to  determine  phase  angle  and  amplitude.   The  theo- 
retical advection  was  determined  by  taking  the  mean  zonal 
component  of  the  wind  (u)  over  a  given  time  period.   A  per- 
centage of  the  actual  advection/theoretical  advection  vs 
latitude  for  a  G-hour  period  is  shown  in  figure  15. 

Experiment  IV.   This  experiment  was  designed  to  analyze 
the  3  effect  as  given  in  equation  (4.2).   It  can  be  shown 
that  the  $  effect  is  represented  by  the  last  term  in  equa- 
tion (4.2).   The  mean  flow,  B  in  equation  (4.2),  was  set  to 
zero  which  eliminates  the  large  N-S  pressure  variation.   The 
8  effect  is  most  pronounced  for  low  wave  numbers.   Wave  num- 
bers 4,  8  and  12  were  run.   According  to  equation  (4.2)  with 
B  set  to  zero,  wave  number  4  should  retrogress  at  24  degrees 
per  day,  wave  number  8  should  retrogress  at  8  degrees  per  . 
day  and  wave  number  12  should  retrogress  at  4  degrees  per 
day.   The  observed  retrogressions  are  shown  in  figures  16, 
17  and  18.   The  actual  retrogressions  were  somewhat  less 
than  expected  due  to  the  mean  divergence  which  reduces  the 
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magnitude  of  g.   The  analyzed  and  48-hour  forecast  surface 
pressure  fields  are  shown  in  charts  A  through  F. 

Experiment  V.   Wave  numbers  4,  8  and  12  were  run  with 
the  mean  flow  adjusted  to  give  a  phase  speed  of  10  degrees 
per  day.   As  the  wave  number  increases  the  second  term  in 
equation  (4.2)  has  less  contribution  and  the  phase  speed  ap- 
proaches the  initialized  phase  speed.   Figures  19,  21  and  23 
show  the  observed  phase  angle  vs  latitude  for  wave  numbers 
4,  8  and  12  respectively.   As  expected,  the  phase  speed  in- 
creased as  wave  number  increased.   The  wave  amplitude  vs 
latitude  are  shown  in  figures  20,  22  and  24.   In  all  three 
cases  the  wave  form  remained  essentially  constant  over  the 
48-hour  forecast  period.   The  analyzed  and  forecast  surface 
pressure  fields  are  shown  in  charts  G  through  L. 
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Figure  11. 


Phase  angle  (degrees  longitude)  vs  latitude,  for 
16-point  N-S  grid,  wave  number  4,  phase  speed 
10°/day  and  A  =  7.0  *  107 .   (Latitudes  with  zero 
wave  amplitude  are  not  included  and  time  is 
given  in  hours.  ) 
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Figure  12. 


Phase  angle  (degrees  longitude)  vs  latitude  for 
16-point  N-S  grid,  wave  number  8,  phase  speed 
10°/day  and  A  =  1.6  *  108 .   (Latitudes  with  zero 
wave  amplitude  are  not  included  and  time  is 
given  in  hours. ) 
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Figure  13.   Phase  angle  (degrees  longitude)  vs  latitude  for 
16-point  N-S  grid,  wave  number  12,  phase  speed 
10°/day  and  A  =  7.0  x  107.   (Latitudes  with  zero 
wave  amplitude  are  not  included  and  time  is 
given  in  hours . ) 
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Percentage  of  actual  advection/theoretical  advection  vs 
latitude  for  wave  number  4  with  A  =  7.0  x  107,  phase  speed 
10°/day,  Coriolis  =  0,  pressure  gradient  =  0,  vertical 
advection  =  0  and  terrain  pressure  and  temperature  constant, 
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Phase  angle  (degrees  longitude)  vs  latitude  for  46-point  N-S 
grid,  wave  number  4,  phase  speed  -24° /day,  B=0  (no  mean 
flow)  and  A  =  7.0  *  107.   (Latitudes  with  zero  amplitude  are 
not  included  and  time  is  given  in  hours.) 


56 


999 


999 


9tfN 


-45  N 


-45*S 


Chart  A. 


-90*S 


Initial  surface  pressure  analysis  for  46-point  N-S  grid,  wave 
number  4,  phase  speed  -24°/day,  B=0  (no  mean  flow)  and 
A  =  7.0  x  io\ 
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Chart  B.   48-hour  surface  pressure  forecast  for  46-point  N-S  grid,  wave 
number  4,  phase  speed  -24° /day,  B=0  (no  mean  flow)  and 
A  =  7.0  x  107. 
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Figure  17.   Phase  angle  (degrees  longitude)  vs  latitude  for  46-point  N-S 
grid,  wave  number  8,  phase  speed  -8°/day,  B=0  (no  mean  flow) 
and  A  =  1.6  x  108.   (Latitudes  with  zero  amplitude  are  not 
included  and  time  is  given  in  hours.) 
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Chart  C.   Initial  surface  pressure  analysis  for  46-point  N-S  grid, 

wave  number  8,  phase  speed  -8° /day,  B=0  (no  mean  flow)  and 
A  =  1.6  x  108. 
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Chart  D.   48-hour  surface  pressure  forecast  for  46-point  N-S  grid, 

wave  number  8,  phase  speed  -8° /day,  B=0  (no  mean  flow)  and 
A  -  1.6  x  108. 
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Figure  18.   Phase  angle  (degrees  longitude)  vs  latitude  for  46-point  N-S 
grid,  wave  number  12,  phase  speed  -4° /day,  B=0  (no  mean  flow) 
and  A  =  7.0  x  10  .   (Latitudes  with  zero  amplitude  are  not 
included  and  time  is  given  in  hours.) 
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Chart  E.   Initial  surface  pressure  analysis  for  46-point  N-S  grid, 

wave  number  12,  phase  speed  -4° /day,  B=0  (no  mean  flow)  and 
A  =  7.0  x  107. 
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Chart  F.   48-hour  surface  pressure  forecast  for  46-point  N-S  grid, 

wave  number  12,  phase  speed  -4° /day,  B=0  (no  mean  flow)  and 
A  =  7.0  x  io7. 
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Figure  19.   Phase  angle  (degrees  longitude)  vs  latitude  for  46-point  N-S 

grid,  wave  number  4,  phase  speed  10° /day  and  A  =  7 . 0  x  10  . 

(Latitudes  with  zero  amplitude  are  not  included  and  time  is 
given  in  hours.) 
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Figure   20.      Terrain  pressure  amplitude  vs   latitude   for   initial   field 

and   48-hour   forecast,    wave  number   4,    phase  speed   10° /day 
and  A  =   7.0  x   107. 
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Chart  G.   Initial  surface  pressure  analysis  for  46-point  N-S  grid,  wave 
number  A,  phase  speed  10° /day  and  A  =  7.0  *  107. 
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Chart  H.   48-hour  surface  pressure  forecast  for  46-point  N-S  grid,  wave 
number  A,  phase  speed  10°/day  and  A  ■  7.0  x  107. 
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Figure  21.   Phase  angle  (degrees  longitude)  vs  latitude  for  46-point  N-S 
grid,  wave  number  8,  phase  speed  10°/day  and  A  =  1.6  *  108. 
(Latitudes  with  zero  amplitude  are  not  included  and  time  is 
given  in  hours.) 
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Figure  22.  Terrain  pressure  amplitude  vs  latitude  for  initial  field 
and  48-hour  forecast,  wave  number  8,  phase  speed  10° /day 
and  A  =  1.6  x  io8. 
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Chart  I.   Initial  surface  pressure  analysis  for  46-point  N-S  grid, 
wave  number  8,  phase  speed  10° /day  and  A  =  1.6  x  10 e. 
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Chart  J.   48-hour  surface  pressure  forecast  for  46-point  N-S  grid, 
wave  number  8,  phase  speed  10° /day  and  A  =  1.6  *  10  . 
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Figure  23.   Phase  angle  (degrees  longitude)  vs  latitude  for  46-point  N-S 

grid,  wave  number  12,  phase  speed  10°/day  and  A  =  7.0  x  107. 

(Latitudes  with  zero  amplitude  are  not  included  and  time  is 
given  in  hours.) 
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Figure  24.   Terrain  pressure  amplitude  vs  latitude  for  initial  field 
and  48-hour  forecast,  wave  number  12,  phase  speed  10° /day 
and  A  =  7.0  x  107. 
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Chart  K.   Initial  surface  pressure  analysis  for  46-point  N-S  grid, 
wave  number  12,  phase  speed  10° /day  and  A  =  7.0  *  10  . 
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Chart  L.   48-hour  surface  pressure  forecast  for  46-point  N-S  grid, 
wave  number  12,  phase  speed  10° /day  and  A  =  7.0  x  107. 
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VI.   CONCLUSION 

The  model  remained  well-behaved  throughout  all  experi- 
ments and  the  phase  speeds  were  less  than  for  the  non-diver- 
gent model,  however,  they  equaled  or  exceeded  theoretical 
expectations  as  described  by  Williams  (1972).   The  model,  as 
it  is  presently  constructed,  appears  to  have  the  desired 
horizontal  flexibility  as  evidenced  by  the  fact  that  the  ex- 
periments were  carried  out  with  different  grid  spacings. 
The  vertical  flexibility,  although  built  in,  was  not  tested 
in  this  report.   The  simplified  nature  of  the  model  pre- 
cluded direct  comparison  with  existing  global  models  such  as 
FNWC!s  5- level  model. 

Future  research  efforts  with  this  model  should  most 
likely  include  expansion  to  5-levels,  treatment  of  cross 
polar  flow  and  initialization  with  real  data. 
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